function [V,Pa,F] = ValueFunction(States,pproc,theta,opt)
% compute value function, action probabilities and unconditional age
% distribution

% Before beginning the estimation per se, we need the matrix with exogenous 
% transition probabilities
Regions = unique(States.RegionsForExoStates);
Nr = length(Regions);
for k = 1:Nr % State transition:
    reg = Regions(k);
    [F{reg}] = StateTrans(States.AllStatesMap,pproc.P{reg});
end

nstates = size(States.AllStatesMap,1);

% Starting value:
V    = ones(nstates,Nr); % One value function vector for each region
Vnew = V;
max_iter = 10000; tol = 10^-4;


for k = 1:Nr
    reg = Regions(k);
    [U{reg},~] = Payoff(States.AllStatesMap,States.MidPointsFixedStates,States.FixedStatesMap,States.SRDMidPoints,theta,pproc.y{reg});
    iter = 0; err = 1;
    while err > tol && iter < max_iter

        iter = iter + 1;

        %             % Outside sugar, there is no third choice:
        %             Vnew_temp  = log(exp(U{r}(:,1) + opt.rho*F{r}{1}*V(:,r)) + exp(U{r}(:,2) + opt.rho*F{r}{2}*V(:,r)));
        %             Vnew(V0,r) = Vnew_temp(V0);

        Vnew(:,k) = log(exp(U{reg}(:,1) + opt.rho*F{reg}{1}*V(:,k)) + exp(U{reg}(:,2) + opt.rho*F{reg}{2}*V(:,k)) + exp(U{reg}(:,3) + opt.rho*F{reg}{3}*V(:,k)));
        % Compute err:
        err = max(abs(Vnew(:,k)-V(:,k)));
        V(:,k) = Vnew(:,k);
        
        if iter == max_iter
            disp(' Max iter reached at Value function iteration ')
        end
    end
    
    [Pa{reg}] = ChoiceProb(F{reg},U{reg},V(:,k),opt);

end

end


